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ABSTRACT 

Scalar fields, strongly coupled to matter, can be present in nature and still be invisible to local 
experiments if they are subject to a screening mechanism. The symmetron is one such mechanism 
which relies on restoration of a spontaneously broken symmetry in regions of high density to shield 
the scalar fifth force. We have investigated structure formation in the symmetron model by using 
N-body simulations and find strong observable signatures in both the linear and nonlinear matter 
power spectrum and on the halo mass function. The mechanism for suppressing the scalar fifth force 
in high density regions is also found to work very well. 



1. INTRODUCTION 

Our current standard model of cosmology, ACDM, has 
been very successful in explaining a large range of obser- 
vations probing a vast range in length scales. We should 
nevertheless be open for the possibility that ACDM is 
just a first order approximation of some more fundamen- 
tal theory. Many theories of high energy physics, like 
string theory and supergravity, p redict light grav it ation- 
ally co upled scalar fields (see e.g. |Binetruy| ( |2006| ; |Linde] 
(2008) and references therein). These scalars may play 
the role of dark energy (quintessence). If these scalar 
fields have non-minimal coupling to matter fields, then 
they could mediate extra forces which are potentially de- 
tectable in local experiments. 

Over the last decades, several laboratory and solar sys- 
tem experiments have tried to detect a sign of such fun- 



damental coupled scalar fields (fAdelberger 2002 ; Hoskins 



et al. 1985 Decca et al. 2007[ |rjertotti et al 
the results so far show no signature of them? 



2DD3J , 
Naively, 

the results of these experiments have ruled out any such 
scalar fields. However, one should bear in mind that a 
coupled scalar field might exist but is undetected just 
because it is either very weakly coupled or very heavy. 

To this da y we know thre e types of theoretical mech- 
anisms (see |Khoury| (2010) for a review) that can ex- 
plain why such light scalars, if they exist, may not be 
visible to experiments performed near the Earth. One 



such class, the chameleon mechanism (iKhoury & Welt- 

" 1)05| IMot-i 



Brax et al.||2004[ |Clifton et alW 
!004a|b ) , operates when the scalar; 



ma EI |20(H _ 

arrow||2004a|pp , operates when the scalars are cou- 
pled to matter in such a way that their effective mass 
depends on the local matter density. In space, where the 
local mass density is low, the scalars would be light and 
deviations from General Relativity would be observed. 
But near the Earth, where experiments are performed, 
the local mass density is high and the scalar field would 
acquire a heavy mass making the interactions short range 
and therefore unobservable. 
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The second mechanism, the V ainshtein mechanism 
Vainshtein||l972| |Deffayet et~al~1 12002| |Arkani-Hamed 



et al.|2003p , operates when the scalar has derivative self 
couplings which become important near matter sources 
such as the Earth. The strong coupling near sources 
essentially cranks up the kinetic terms, which translates 
into a weakened matter coupling. Thus the scalar screens 
itself and becomes invisible to experiments. This mech- 
anism is central to the phenomenological viability of 



theories ( 


Dvali et al.||2000t |de Rham et ah] 


2008: Nicolis 


et al.[2008; 

Gabadadze 




Hinterbichler et al.||2010l |Mota et al.||2010l 

20091 Ide Rham||2010| IBrax et al.||2011a|). 
I i r 



Fiie last mechanism, the one explored m this paper, 
is the symmetron mechanism ( [Hinterbichler & Khoury 
20l0l IHinterbichler et aT[|201l| ; [Olive k Pospelov||iMOg[ 
Brax et al.||2011b[ ). In this mechanism, the vacuum ex- 
pectation value (VEV) of the scalar depends on the lo- 
cal mass density, becoming large in regions of low mass 
density, and small in regions of high mass density. By 
taking the coupling of the scalar to matter to be propor- 
tional to the VEV, we can have a viable theory where 
the scalar couples with gravitational strength in regions 
of low density, but is decoupled and screened in regions 
of high density. This is achieved through the interplay of 
a symmetry breaking potential and a universal quadratic 
coupling to matter. 

In vacuum, the scalar acquires a VEV which sponta- 
neously breaks the Z2 symmetry (j) — > —(f). In the regions 
of sufficiently high matter density, the field is confined 
near <j> — 0, and the symmetry is restored. The fifth 
force arising from the matter coupling is proportional to 
<f> making the effects of the scalar small in high density 
regions. 

As op posed to Chameleons, where the strongest con- 
straints dMota fc Shaw||2007| |Brax et al.|[2^10bl |2008 
Gannouji et al.|2010||Gies et al.|2008[|Brax et al.|20(J7a|b 



Mota & Shaw|l2006) comes from laboratory experiments 



which in effect, washes out any observable effects in the 
solar system, the symmetron predicts a host of obser- 
vational signatures in experiments designed to look for 
deviations from GR, which are just below the currents 
bounds and within reach of the next generation experi- 



2 



merits. 

The cosmology of coupled scalar field models are usu- 
ally strongly constrained by local gravity experiments, 
which could put limits on the range and the coupling 
strength of the scalar field. There do exist several cases 
in which signatures on the linear perturbations are found, 
but in most cases the range of the field is well below linear 
scales. To proceed into the region of nonlinear structure 
formation one can use the spherical collapse model to ob- 
tain the qualitative behavior, but in order to obtain accu- 
rate quantitative results deep into the nonlinear regime 
one is almost required to perform N-body simulations. 

Studies of coupled scalar field models by using N - body 
simulations (iZhao et al.||2011a| iBrax et al.||2011c| |Li & 



Barrow 2011a 
Zhao 2010 Ba: 



Li et ai.||Mlal|Zhao et aiPTHllLTT 

i et al.||2010: Baldi 2009) have revealed" 



several interesting signatures which can in principle be 
de tected by observation s in the near future. For example, 



Zhao et al. (2011b I it was found that f(R) theories 



can give rise to a dependence on the environment of the 
dynamical to lensing mass ratio of halos; an observable 
feature that is not found in ACDM. 

In this article we will study the effects a symmetron 
field has on structure formation. By performing high 
resolution N-body simulations we demonstrate explicitly 
how the Symmetron mechanism works in screening the 
fifth force and obtain observables as the matter power 
spectrum and the mass function. 

2. THE SYMMETRON MODEL 

The action governing the dynamics of the symmetron 
model is given by 



dx 4 



2 pl 



v{4>) 



(1) 



where g is the determinant of the metric g^ Ul R is the 
Ricci scalar, ^ are the different matter fields and M p \ 
, where G is the bare gravitational constant. The 

matter fields couple to the Jordan frame metric g^ v via 
a conformal rescaling of the Einstein frame metric g au 
given by 

g„„ = A 2 {4>)g^ (2) 

The coupling function A(<f>) is chosen to be an even poly- 
nomial in <p (to be compatible with the <fi — > — <f> symme- 

*«- i +K£) ,+o (£) < 3 > 

described by a single mass scale M. For the range of 
parameters we are interested in we have (j^j 1, 

thus, we can neglect the higher order correction terms. 
The potential is chosen to be of the symmetry breaking 
form 

1 _ _ 1 

(4) 



V(<f>) = V - \p 2 <t> 2 + \M 4 



where Vq is a cosmological constant (CC). We will for 
simplicity absorb all contributions to the CC into Vo by 
simply putting Vq = A. We will later see that A must be 
taken to be the usual CC to obtain late time acceleration 



of the Universe. The field equation for (f> follows from the 
variation of the action Eq. (Ill with respect to (j> and reads 

□0 = VeBJ (5) 

The effective potential is given in terms of the trace, T m , 
of the matter energy-momentum tensor by 



T 



1 / Pm 

2 \M 2 



f U 2 + W 



2 + 4 A ^ 



(6) 
(7) 



depending on the two mass scales //, M and the dimen- 
sionless coupling constant A. It is convenient to define 
the critical matter density (and the critical redshift) 



Pssb = n 2 M 2 = 3iJ 2 M^ m (l + z SSB ) 



(8) 



where VL m is the matter density parameter in the Uni- 
verse today and Hq is the Hubble parameter. In regions 
where p m > pssb (where p m is the local matter density) 
the symmetry —>•—(/> is upheld and the effective po- 
tential has a minimum at m i n = 0, whereas in regions 
where p rn < psSB the symmetry is spontaneously broken 
and the field acquires a VEV 



= ±<Wi- 



PSSB 



(9) 



where (f>Q = is the symmetry breaking VEV for p m — > 

0. The mass of small fluctuations around the minimum 
of the effective potential is given by 

( P 



ml = I4ff . 



- 1 ) p z + 3\(fj 

V PSSB / 
^-l), Pm>PSSB 
1 _ ^) ' ^ < PSSB 



2 

min 



(10) 



The symmetron field acquires the longest range, 



in low density regions where 



An 



1 



(11) 



For future convenience we introduce the dimensionless 
quantity L = M p°/ h , which is the maximum range of the 

symmetron mediated force in units of Mpc//i. 

The gravitational field equation for g^ v is given 
by 



Guv — SirGTn^ 



(12) 



where the total energy-momentum tensor T^u is the sum 
of the matter and scalar field parts: 



A{4>)TZ + - g uv -{d<P) z + V(4>) (13) 



1 



Note that the matter part itself is not conserved, but 
instead satisfies 



rrv,. 



(14) 
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In N-body simulations we are interested in describing 
the matter sector by particles and the energy-momentum 
tensor of an individual particle with mass mo at position 
ro is given by 

T^(r) = ^=S(r-r )r^ (15) 



-fl 

where r is the ge ner al spatial coordinate. Taking the 
divergence of Eq. ( 12 I and using the Bianchi identity we 
get the geodesic equation for the matter particles 



' ' 1 a~{' ' — 



dlogA(^) 



(16) 



which for A = 1 reduces to the standard geodesic equa- 
tion in general relativity. 

2.1. The Symmetron Mechanism: Local Constraints 



From Eq. ( 16 1 we see that the symmetron field gives 
rise to a fifth iorce on the matter fields which, in the 
nonrelativistic limit, is given by 



M 2 



V<t> 



P ( <$> 



M, 



pi 



V0 



(17) 



where we have introduced the coupling constant /3 = 

0oM p i 
M 2 ■ 

The static spherical symmetric solutions of the field 
equations were found in Hinterbichler & Khoury (20101. 
For two test masses in a region where <p = <pB it was 
shown that that the fifth force is simply 



= 2/3 2 



(18) 



In a low density region (p <C Pssb) we have 4>b = (f>o & n d 
the fifth force is comparable with gravity for j3 = 0(1). 

For very large bodies, the situation is quite different. 
The symmetry is restored in the interior of the body and 
the fifth force on a test mass outside becomes 



F N 



2f3 2 



a ' 



At, 



PSSB 



Pbody \ -^body 



(19) 

The fifth force is suppressed by a factor a -1 -C 1 — sim- 
ilar to the thin shell factor found in chameleon theories 
( |Khoury fc Weltman|[2004 l. 

We also see that if the test masses are inside a screened 
region <C 1^ the force will be further suppressed. 

Since the field is long ranged (and universally coupled) 
in almost all situations today the theory is best con- 
strained by solar system experiments which have been 
performed with high precision. 

It turns out that as long as our Galaxy is sufficiently 
screened (10 < ac), our Sun will also be screened and 
the combined effects discussed above are enough to evade 
the current parameterized post Newtonian (PPN) con- 
straints. 

By assuming that 4> — > <po outside our Galaxy, i.e. that 
our galactic neighborh ood is not screened, these con- 



straints were derived in Hinterbichler & Khoury (20101; 



Brax et al. (2011b I, and require 
M < 10~ 3 



pi 



(20) 



If the assumption about the value of <fi outside our 
Galaxy, which is very likely to be true, can be relaxed 
then the bound above can be relaxed somewhat as well. 
The constraint on M turns into a constraint on the range 
of the field and the redshift in which the SSB takes place: 



A < 2.3 




(1 + z SSB )- 3/2 Mpc// l (21) 



Thus for transitions that take place close to the present, 
the fifth force can have a range of at most a few Mpc/h. 

2.2. Physical Parameters 

In the rest of this article, instead of working with the 
parameters {p,M, A}, we will instead choose to work 
with the more physically intuitive quantities {L, f3, zssb} 
: the cosmological range of the fifth force in Mpc//i, the 
strength of the fifth force relative to gravity and the red- 
shift at which the SSB takes place in the cosmological 
background. 

The transformation between the two sets of parameters 
is given by 



I 1 


2998 


Ho 


V2L 


M 






io- 3 v 




A = 


/10 60 





(L27 



H 



pi 



0.27 



L 
2^36 

2 



(1 + ^ssb) 



3/2 



1.38-10 



-100 



/32L6 (1 + Zssb)C 



(22) 



(23) 



(24) 



For typical parameters L ~ f3 ~ 1 and zssb ~ we have 
fi - 10- 57 M pl , M - 10~ 3 M pl and A - lO" 100 : thus the 
symmetron is very weakly coupled. 

We will choose to work with values of the parameters, 
that are close to the local constraints, and in which the 
symmetron can produce observable cosmological effects. 
This means we will be most interested in the parameter 
space L = 0(1), j3 = 0(1) and < z SSB < 2. 

3. SYMMETRON COSMOLOGY 

In this section we discuss the cosmological evolution of 
the symmetron field from the background evolution to 
linear perturbations and derive the nonrelativistic limits 
of the field equations to be implemented in the N-body 
code. The analysis in this section is mainly for com- 
parison with the N-body simulations. For a more thor- 
ough discussion regarding the background c osmology and 
linear pertu rbati ons in the symmetr on see |Hinterb ichlcr 
et al. (20111 and Brax et al. (2011b I respectively. 



3.1. Background Cosmology 

The background evolution of the symmetron in a flat 
Friedmann-Lemaitre-Robertson- Walker (FRLW) metric 

ds 2 = -dt 2 + a 2 (t)(dx 2 + dy 2 + dz 2 ) (25) 
is determined by the field equation 

4> + 3Hj> + V cfi j = (26) 
together with the Friedman equations 



3ff 2 M p 2 ! = Pm A(4>) + , 

p rn + 3Hp m = 



(27) 
(28) 
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Fig. 1. — The background evolution of the symmetron for /3 = 1, 
L = 1 and zssb = 2 together with the analytical background 
(dashed lines). The symmetry is broken at z = 2 and the field 
settles at one of the two branches. 



where 



= A- 



■ 2 ±2 



A</> 4 



(29) 



When the field follows the minimum of the effective po- 
tential we have 



A 



A 



< 



= /3 



2 PSSB 

A 



~ AA 

< 10- 6 /? 2 (l + Z S SB 



M 



(30) 
(31) 



thus for f), zssb ~ 0(1) the dynamical part of the poten- 
tial is too small to contribute significantly to the energy 
density of the Universe and we are left with the cosmo- 
logical constant to account for dark energy. 
In the same regime, the coupling function A(<j>) satisfies 



\A(4>) 1| 



Al 



</3 2 



M 



< l(T 6 /3 2 (32) 



which is also too small to produce an observable effect on 
the background expansion. This implies that the sym- 
metron evades Big Bang Nucleosyntesis (BBN) bounds 
on the variati on of m asses of the standard model parti- 
cles (see Sec. (3.5.11). It might be possible to make the 
symmetron responsible for dark energy by changing the 
form of the potentia l and coupling. One such modifica- 
tion was proposed in |Hinterbichler et al. ( 2011[ ), however 
it was shown that tuning of the parameters was required 
to yield the desired late time cosmology. 

In Fig. ([lj we see the background evolution <p(z) for 
z ssb = 2 together with the analytical minimum. No- 
tice that the field does not immediately start to follow 
the minimum right after SSB. This has important conse- 
quences for the evolution of the perturbations which will 
be discussed in section l3~4l 



3.2. Linear perturbations 

The most general metric in a perturbed FLRW space- 
time is given by 

ds 2 = -(1 + 2a)dt 2 - 2aB, l dtdx i (33) 

+a 2 ((1 + 2ip)5ij + 27,iy) dx i dx j 

where the covariant derivative is given in terms of the 
spatial metric, which in the case of our flat background 



reduces to 6ij. We decompose the field cj> into the back- 
ground and perturbated parts: (j)(x,t) — cj>(t) + 5cf>(x,t). 
The energy-momenntum tensor of nonrelativistic matter 
can be decomposed as 



rpO _ 

1 o — 



-Pi 



l (l + <5 m ), T? = - Pn 



V,i 



(34) 



where v is the peculiar velocity of nonrelativistic matter 
and S m is the matter density perturbation defined by 



S m = 



5p„ 



Pn 
Pn 



Sp r , 
Pm 



in the comoving gauge (35) 



The equation determining the evolution of the pertur- 
bations, neglecting anisotropic stresses, follows from the 
Einstein equations. The scal ar perturbations can be read 
off from the formulation of (|Hwang fc Noh| |2002[ ) , which 
is independent of gauge. In the following we use units of 
M p i = 1. After solving for the different metric potentials 
we find that the scalar perturbations, in the comoving 
gauge {v — 0) are determined by 



2HS„ 



(36) 



AP 



6H 2 + 6H + fl m H 2 - 



'2 



2(j> 



AP 



V eSA ( 1 + — 2 ) 60 - 2<f>5<t> = 




(37) 



In studying the perturbations it is convenient to intro- 
duce the growth index 



7(z,fc) 



log(fi m (z)) 



(38) 



In ACDM we have 7 « 0.55 (for 0.2 < Q m < 0.3), 
which is scale and almost redshift independent. In 
modified theories however, 7 can have significant scale 
and r e dshift dependence as s hown in |Gannou,ji et aL| 



|2010t;|Tsujikawaet al.l ( |2009[) ; 
Mota & Wintherl a2011DTorTh. 



ota fc Winther| ( |201l| | 
models. 



Brax et al ,| ( |2008[|2010'aT - 

e case of chameleon// (ft) 



If we assume that the field is rolling slowly along 
the minimum we can neglect all terms proportional to 
<p and the oscillating term V c q^. The perturbations in 
4> will evolve more slowly than the perturbations in 5 m 
for scales deep inside the Hubble radius, thus, the term 

Pm(3,<f,S m and (m^ + ^i)5cf) will dominate over the S(f> 
time derivatives in Eq. (37 1. Under these assumptions, 



5 



we can simplify Eq. ( 36 I to 
5" 



2H8 m 
Gcs 



= i 



3 o n 2 GeS a 



(39) 
(40) 



which are the equations we use to integrate the pertur- 
bations. 

At the time before SSB we have <j> w and therefore 
= 1. After SSB the field approaches the minimum 



«o = m this regime we have 



~G~ 



2(3 



2 



(41) 



Thus small scales will feel a stronger gravitational con- 
stant. 

In Fig. (pi we show the redshift evolution of 7 for 
several different wave numbers and in Fig. Q we show 
contour plots for 7(2 = 0) for two comoving wavenum- 
bers. 

The growth rate on really large scales (k < 0.01/i/Mpc) 
is not affected by the symmetron fifth force unless 
L,f3 ^> 1. However on the smallest, linear scales we 
can still have a deviation from the predications of 
GR. Note that we have integrated the perturbations 
using the approximation Eq. (39) instead of the full 
equatio ns E q. ( 36|37 ). The explanation for this is given 
in Sec. JOb. 



3.3. Linear power spectrum and the CMB 

In Fig. Q we show the factional difference of the lin- 
ear power spectrum of the symmetron to that of ACDM, 



defined as 



AP(fc) 
p(k) 



P(k)-P AC OM(k) Notice that on 

r\CDM (K) 



ear scales (k < 0.1/i/Mpc) the power spectrum is very 
close to ACDM. Going down to scales comparable to the 
length scale of the symmetron (fc ~ j-h/Mpc) the power 
spectrum starts to deviate significantly. However, in this 
regime the perturbations are already nonlinear and we 
cannot trust the results of the linear perturbation theory. 
Once we discuss the N-body results we will see that the 
symmetron mechanism is at work in this regime, thereby 
suppressing the predication of linear perturbation theory. 

The relative short range of the fifth force means that 
it will not affect the CMB unless L,/3 > 1. Take L = 1 
and f3 — 2 as an example: we find a maximal increase in 
power (due to the ISW effect) of ~ 0.25% for multipoles 
around I ~ 100. One needs a much larger j3 and/or L to 
have a detectable signature in the CMB. The second case 
is not allowed by local experiments while the first case 
implies a growth rate of the linear perturbations which 
should have difficulty satisfying constraints coming from 
large scale structure surveys. 

A more thorough analysis of the linear p erturbations 
in the symmetron model can be found in |Brax et al.| 
(2011b I. There it was shown that strong signatures ap- 



pears m other interesting linear observables such as the 
weak lensing slip parameter and the modified gravity pa- 
rameter. 

3.4. Tachyonic Instability in the perturbations 



The perturbations in Sec. ( |3.2|) w ere integrated using 
the approximate equation Eq. (39 I, which is equivalent 



to using the analytical minimum as the background field. 
The reason we did not use the full equations is because 
perturbations theory breaks down close to zssb- 

Immediately before zssb the field is still close to <\> = 0; 
as z — > zssb the mass of the field vanishes. This means 
that the field cannot follow the minimum and starts to 
lag behind as seen in Fig. (111. The global minimum of 
the effective potential (f> = Oiiow becomes a local maxi- 
mum and the mass squared of the field becomes negative 
leading to an exponential growth in the pert urb ations. 

To see what happens we can simplify Eq. (37 1 by dis- 
carding all but the most important terms. 



+ m J + — ~ 



If 



(42) 



< then the solution to the above equation 



reads 8(f) oc e*v / ' m '* +fe2 / a2 ' which is exponentially growing. 

In a realistic situation the field would roll very quickly 
down from the false minimum <fi = 0, making posi- 
tive, and thereby stabili zing the field close to the symme- 
try breaking minimum ( |Felder et al.|2001[ ). Since pertur- 
bation theory is only valid as long as the perturbations 
8(f) are small we get a breakdown of the perturbation the- 
ory when using the true background solution. The blow 
up in 8(f), in turn, leads to a blow up in the matter per- 
turbations and the numerical results cannot be trusted. 

We have investigated this further by using N-body 
simulations. In Fig. (10 1 we see a snapshot of the <t>- 



distribution both before and after z — zssb = 2.0. There 
we see the same sort of behavior as is familiar from sym- 
metry breaking in condensed matter physics: symmetry 
breaking takes place at different places at different times 
according to the local matter density. This type of dy- 
namics is not taken care of in the standard perturbation 
theory approach which leads to the apparent instability. 

Note that by using the analytical minimum when in- 
tegrating the perturbations we do not have control over 
the accuracy of our results. A full analysis of this phe- 
nomenon could be handled with N-body simulations, but 
in our simulations we have not explicitly taken into ac- 
count the time variation of the scalar field (we work in 
the quasi-static limit) and our simulation box is also too 
small to reach far enough into the linear regime. We 
leave this study for future work. 

3.5. Varying Constants 

3.5.1. WMAP constraints on particle mass variation 

One important constraint on coupled scalar field 
theories comes from time variation in the gravitational 
constant G in the Jordan frame, or equivalently in the 
masses of the standard model particles in the Einstein 
Frame. WMAP constrains any such variation to be les s 
than about 5% since recombination ( |Nagata et al.[ 2004 ) . 
Light element abundances provide similar co nstraints 
between th e time of nucleosynthesis and today ( |Accetta| 
eTaT1p90 l. 



Due to the conformal coupling to matter, A((f)), a 
constant mass scale in the Jordan frame becomes time 
and space varying in the Einstein frame. The mass 








Fig. 3. — The growth index f(z = 0) for 2sSB = 1-0 (left) and 2ssb = 2.0 (right) for two comoving wavnumbers: k = 0.2Mpc/h (above) 
and k = O.OlMpc/h (below). The red region shows the GR regime 7 ~ 0.555, the blue region show the regime where 0.5 < 7 < 0.55, the 
green region shows 0.4 < 7 < 0.5 and the white region shows 7 < 0.4. 
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Fig. 4. — The linear power spectrum relative to that of ACDM 
for three different SSB redshifts: zgsB = 0.5 (dotted), zssb = 1-0 
(dashed) and zssb = 2.0 (solid). We have fixed L = I and shown 
the results for the two values (3 = 0.5 and /3 = 1.0. 

variation between today and recombination is given by 

Am _ A((j) ICC ) - A(<^today) ^ 1 ( 0today\ 2 

~ 2 V M ) ^ ' 



A{(f>rec) 



where we have put <p rcc ~ since z rcc 3> ^ssb in all 
interesting cases. If we further assume c^today = 4>o we 
get the conservative constraint 



6 fl 2 



(44) 



Am^l =B 2 (M_ 

m 2\MJ P \M pX/ 

The WMAP constraint | ^ | < 0.05 is satisfied for all 
/3<100. 

3.5.2. Fine structure constant 

Analysis of absorption spectra of quasars have led some 
to claim that the fine structure constant a might have 
evolved by approximately one part in 10 5 over the red- 
shift range 0.2 < z < 3.7. If this turns out to be true, 
then general covariance would imply that a can vary both 
in space and in time, that is, it must be a function of a 
field. 

Since we have so far assumed that the symmetron cou- 
ples conformally to matter fields, and since the Maxwell 
action is conformally invariant, at tree level the sym- 
metron does not lead to a time varying a. By considering 
a coupling of the symmetron to photons of the form 



where 



Ayfa) = 1 + 



M 



(45) 



(46) 



then variations in (f> will lead to variations in a. Here 
£ 7 is the symmetron-photon coupling relative to the 
symmetron-matter coupling. The variation in the fine 
structure constant between Earth (E) and another place 
(S) in the Universe is given by 



Aq 

a 



A^e) ~ 2 



<Ps_ 
M 



(47) 



If S is a very low density environment where <ps 
then 



a 



C 7 io 



-6^2 



(48) 



which for 0(1) < j3, C 7 is close to the reported detection. 

However, the local density in most Ly-a emitting sys- 
tems is usually much larger t han the cosm ological back- 
ground density today (see e.g. Brax et al. ( |2004| and ref- 
erences therein) , which implies <ps <C <po and the above 
estimate becomes even smaller. 

To be able to account for the reported claims we need 
^SSB to be well before the observed redshift of these sys- 
tems and/or these systems to be located in voids to pro- 
duce the desired 10~ 5 effect. This makes it possible that 
the symmetron is responsible for the claimed variations, 
but most likely it will require a fin e tuning C-v ^ 1 - A 
more detailed analysis, as done in Li et al. (2011b I, is 
required to see if this is the case. This is beyond the 
scope of this paper. 



To 



3.6. N-body Equations 
implement the general relativistic equations 



Eq. ( 5|12|13|16 1 in N-body simulations, it suffices to work 
in the nonrelativistic limits, since the simulations only 
probe the weak gravity regime and small volumes com- 
pared with the cosmos. We write the perturbed metric 
in the (flat) conformal Newtonian gauge as 



ds 2 = -a 2 (l + 2E)dr 2 + a 2 (l - 2^)dx l *dx f , 



(49) 



where r is the conformal time and 
coordinate. 



is the comoving 



The scalar field equation of motion in terms of the 
perturbed quantities becomes 



-(1-2H) 



■ - <j)' (2H(l - 2E) — 3' — 3*') 



( RlIL. 

Am 2 



(50) 



Taking the quasi-static limit of this equation, in which 
we can neglect terms such as 3', ty' and H<fi' since the 
time derivative of a quantity is much smaller than its 
spatial gradient, and removing the background part we 
obtain 



M 2 pl 



(p m (f> - p m 4>)(j>) 



M 2 pl 

where we have also used the approximation A(cf>) 
simplify the equation further. 



(51) 
1 to 



The (0, 0)-component of the Ricci tensor and the 
trace of the total energy-momentum tensor in the 
perturbed quantities becomes 



a 2 R° n w -V 2 3 + 3 



- H 2 (1 - 25) 



- 3*" - 3H(5' + 

T « -A{4>)p m - iV(<f>) + 2*) 



(52) 
(53) 
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The (0, 0)-component of the Einstein equation with the 
background part removed gives the nonrelativistic Pois- 
son equation 



V 2 ^ a AttG (p m - p m ) a 3 



(54) 



where we have neglected the contribution from the 
potential (V(<p) — V(cf>)), put A{4>) ps 1 and taken 
$ = aS for convenience. 

The equation of motion for the N-body particles 
follows from the geodesic equation and reads 



2iix 



1 



a 2 M 2 



M 2 ' 



(55) 



By rewriting this equation in terms of the conjugate mo- 
mentum to x, p = a 2 x, we have 



P 

dx p 

~dt ~ 
dp 

~dt 



a- 



1 <b , 



(56) 
(57) 



The equations Eq. ( 51|54|56 l are all we need to put into 
the N-body simulation code in order to study structure 
formation in the nonlinear regime. 

4. N-BODY SIMULATIONS 

Below we describe the algorithm and model specifica- 
tions of the N-body simulations we have performed. We 
also give results from tests of the code to show that the 
scalar field solver works accurately. 

4.1. Outline 

For our simulations we have used a modified version 



of the publicly available N-body code MLAPM (Knebc 
et al. 20011. The modi fications we have ma de follow the 
detailed prescription of |Li fc Barrow| ( 2011a[ ), and here we 
only give a brief description. The MLAPM code has two 
sets of meshes: the first includes a series of increasingly 
refined regular meshes covering the whole cubic simu- 
lation box, with respectively 4,8,16,. ..,Ad cells on each 
side, where Nd is the size of the domain grid, which is the 
most refined of these regular meshes. This set of meshes 
are needed to solve the Poisson equation using multigrid 
method or fast Fourier transform (for the latter only the 
domain grid is necessary). When the particle density in 
a cell exceeds a pre-defined threshold, the cell is further 
refined into eight equally sized cubic cells; the refinement 
is done on a cell by cell basis and the resulting refinement 
could have arbitrary shape which matches the true equal 
density contours of the matter distribution. This second 
set of meshes are used to solve the Poisson equation using 
the linear Gauss-Seidel relaxation scheme. 

The symmetron field is the most important ingredient 
in the model studied here, and we have to solve for it to 
obtain detailed information about the fifth force. In our 
N-body code, we have added a new scalar field solver. 
It uses a nonlinear Gauss-Seidel scheme for the relax- 
ation iteration and the same criterion for convergence 
as the default Poisson solver in MLAPM. But it uses V- 
cycle instead of the self-adaptive scheme in arranging the 
Gauss-Seidel iterations. 



The modified Poisson equation is then solved using 
nonlinear Gauss-Seidel relaxation on both the domain 
grid and the refinements. With the gravitational poten- 
tial $ and the scalar field 4> at hand, we can evaluate 
the total force on the particles and update their mo- 
menta/velocities which are used to advance the particles 
in space. 

4.2. Simulation Details 

The physical parameters we use in the simulations are 
as follows: the present dark-energy fractional energy den- 
sity Q A = 0.733 and Cl m = 0.267, H a = 71.9km/s/Mpc, 
n s = 0.963 and a% — 0.801. We use a simulation box 
with size 64Mpc//i, in which h = iJ /(100km/s/Mpc). 
We simulate 9 different models, see Table, (flj) for the 
symmetron parameter values. 

These parameters are chosen so that they predict local 
fifth forces which are of the same order of magnitude as 
allowed by current experiments and observations and are 
such that we can see the effect of the different parameters. 
Note that the energy density in the symmetron is always 
much less than that of dark energy and therefore does not 
alter the background cosmology which in all runs will be 
that of ACDM. 

In all those simulations, the particle number is 256 3 , 
so that the mass resolution is 1.114 x 10 9 Mpc//i. The 
domain grid is a 128 x 128 x 128 cubic and the finest re- 
fined grids have 16384 cells on each side, corresponding 
to a force resolution of about 12kpc//i. The force resolu- 
tion determines the smallest scale on which the numerical 
results are reliable. Our simulations are purely N-body, 
which means that baryonic physics has not been included 
in the numerical code. 

4.3. Initial conditions 

Initial condit ions for the simulation was generated us- 
ing GRAFIC2 ( |Bertschinger|2001[|Pruriet et al.|200"8| by 
using the parameters described above, but where we also 
included baryons with a density parameter fib = 0.045 
(and a dark matter density il m = 0.267 — f2&). We use 
the same initial conditions for all the simulations in order 
to see clearly the effect of the symmetron compared with 
ACDM. 

This choice needs some justification. First of all, we 
start the simulation at z = 49, a time in which the sym- 
metron has no effect on the growth of the perturbations. 
This means that the only change the symmetron will 
have on the initial conditions is on the value of cr 8 today 
which is used to normalize the perturbations. Since the 
symmetron field has a rather short range compared to 
the linear regime we do not expect a large effect on erg 
for the range L < 0(1) we are considering. 

To check this assumption we integrated the perturba- 
tions and calculated the value of o~g, (by normalizing to 
the CMB) for our simulation models and found that the 
model with L — 1 that is furthest away from ACDM, 
namely F in which zssb = 2.0, L = 1 and [3 = 1, only 
has <jg 
conditions 

If one is to consider models in which L is much larger 
than 1 then this becomes an issue that should be dealt 
with properly. 

4.4. Code tests 



1.01cr 8 jCDM justifying the use of ACDM initial 
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Model: 


A 


B 


c 


D 


E 


F 


G 


H 


ACDM 


2SSB 


0.5 


0.5 


1.0 


1.0 


2.0 


2.0 


1.0 


1.0 


0.0 




0.5 


1.0 


0.5 


1.0 


0.5 


1.0 


0.5 


1.0 


0.0 


L 


1.0 


1.0 


1.0 


1.0 


1.0 


1.0 


2.0 


2.0 


0.0 



TABLE 1 

The symmetron parameters used in our simulations. 



Before we run simulations we have to make sure that 
the scalar field solver, which is the main modification to 
the MLAPM code, works accurately by performing code 
tests for situations where the outcome is known from 
analytical solutions. 

The scalar field solver uses the nonlinear Newton- 
Gauss-Seidel relaxation scheme to compute x = ^~ > an d 
an indicator that it works is to show that, given the 
initial guess of the solution that is very different from 
the true solution, the relaxation could produce the latter 
within a reasonable number of iterations. We consider a 
simulation box with homogeneous density (obtained by 
putting particles on a regular grid inside the simulation 
box), then the true solution is given by \ — x'- the back- 
ground solution. We therefore make an initial guess for \ 
which is randomly scattered around x an d let the scalar 
field solver solve for x- The results for |% — x\ before and 
after the relaxation scheme are shown in Fig. |5|. The 
differences between the initial guess and the true solution 
varies between 0.001 and 0.1 while after the relaxation 
the difference is of order 10~ 8 . By using double preci- 
sion numbers in all the calculations we obtained exactly 
the analytical solution (to double precision « 10~ 15 ), 
while for using only floating point numbers the accuracy 
dropped to 10~ 6 which is exactly the accuracy in floats. 
This shows that the scalar field solver works accurately. 

The most important effect of the symmetron is the 
screening mechanism in which the local value of the field 
should be pushed down towards \ — in high density 
environments. We therefore consider a spherical over- 
density, located at the center of the box, with a given 
radius R, homogeneous density p c inside R and embed- 
ded in a background of homogenous density pj,. The 
analytical solution reads 



X(r) =x(0) 



X{r) 
where 



~-Xb 



sinh [m c r] 

MR) 



Xb)R 



-m b (r-R) 



r < R 



r>R 



(58) 
(59) 



(m 2 ) ' 




-d) 

1 + rribR 



sinh{m c R) + m bR / 



(60) 



For the trial solution on the grid we use the background 
value Xb and we perform the test for a range of densi- 
ties p c . The results after relaxation for the most mas- 
sive cases are shown in Fig. (Rjl. There are some small 
discrepancy from the analytical solution in the region 
R < r < 2R for the most extreme cases p c > 10 3 ph. This 



i 




x (Mpc/h) 

Fig. 5. — The scalar field relative to the analytical solution before 
(random initial values above) and after (below) the Newton-Gauss- 
Seidel relaxation. 




r/R (R= Mpc/h) 

Fig. 6. — The scalar field value as function as distance from 
the center for a spherical over-density embedded in a background 
of homogenous density pi, together with the analytical solution 
for p c = 4000, 6000 and 8000 times pi,. The points shown here 
are calculated by binning the scalar field value using a bin-width 
A(r/i?) = 0.01 and taking the average. We used the same amount 
of particles, 128 3 , in each run so that the background density p;, 
differs for the three cases shown above. 

is not a surprise as the density suddenly drops over 3 or- 
ders of magnitude at r = R, meaning that we need a lot 
of particles in this region in order to get accurate results. 
In the region r < R and r > 2R the scalar field solver 
produces the analytical solution to high accuracy. 

5. NUMERICAL RESULTS 

In this section we present the results from the simu- 
lations, including the snapshots, the matter power spec- 
trum and the halo mass function. 

5.1. Snapshots 
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In the symmetron model \ = ^ > an d thereby the fifth 
force, is suppressed in high density regions. In this sub- 
section we demonstrate these qualitative features using 
some snapshots. 

Fig. ([7 1 shows the ratio of the fifth force to gravity 
today for redshift both before and after zssb- 

At early times, the density is high everywhere and we 
expect the fifth force on all particles to be strongly sup- 
pressed. At later times we expect a screening in regions of 
high matter density. These predictions are confirmed in 
Fig. 0. We see that fifth force on the particles which feel 
a strong gravitational force (i.e. particles in a high den- 
sity environment) is highly suppressed whereas the fifth 
force on particles which feel a weak gravitational force 
(i.e. particles in a low density environment) follows the 
unscreened theoretical prediction ~ 2f3 2 xl(z)F glav i ty 
(see Eq. (f39l) 



Fig. 



Ol shows the density and scalar field dis- 



tribution in a slice of the simulation box at different 
redshifts for the three cases zssb = 0.5, 1.0 and 2.0 with 
/3 = L = 1.0 fixed. 

For redshifts z > zssb, X is very close to the min- 
imum x = almost everywhere in space except in voids 
where the symmetry has already been (weakly) broken. 
When we go down to redshifts z < zgsB the symmetry 
is broken in most parts of the box, except in the high 
density regions where we still have x ~ 0- Comparing 
the scalar field distribution today for runs with different 
^ssb, we see that the earlier the symmetry breaking 
takes place the larger the part of the box which is 
unscreened (x ~ 1) today becomes. This is because 
the critical density for the symmetry breaking is larger 
for larger zssb and therefore the halos have to be more 
massive in order to be effectively screened. 

5.2. Matter power spectrum 

The nonlinear matter power spectrum is an important 
observable and could be used to distinguish among differ- 
ent models of structure formation. As we have seen above 
the symmetron can have a strong effect on the growth 
rate of the linear perturbations for parameters that are 
allowed by local experiments. We expect these signatures 
to sho w up in the nonlinear matter power spectrum. 

Fig. ( 12 1 displays the fractional difference in the mat- 



ter power spectrum from that of ACDM, defined as 
(P(k) - PACDAi(k))/P(k), and in Fig. (ITT} we show the 
actual power spectrums for the symmetron and ACDM 
together with the corresponding predictions from linear 
perturbations theory. 

The power spectrum agrees with the predictions of lin- 
ear perturbation theory on large scales (k < 0.1 Mpc//i), 
but on smaller scales the results found here are weaker 
than the prediction of linear perturbation theory seen in 
Fig. Q. This is because when linearizing the field equa- 
tion we are basically using the background matter den- 
sity everywhere and therefore preventing the symmetron 
mechanism to take effect in suppressing the fifth force 
when matter perturbations become large. In contrast, 
the N-body simulation avoids this approximation by tak- 
ing full account of the suppression of the fifth force. 

The fractional difference relative to ACDM is growing 
with zssb and (3 as the fifth force has more time to 



operate and is stronger. Comparing runs with the same 
P we see an important effect if the symmetry breaking 
is earlier. When zssb = 2.0 the fractional power is 
increasing until we reach a scale where the screening 
mechanism becomes stronger and then starts to decrease 
again towards ACDM, only to start growing again at 
even smaller scales. This is because the critical density 
for having screening is much higher for larger zssb so 
that most halos (which are on small scales and of low 
mass) are unscreened. 



In Fig. ( 13 1 we show the redshift evolution of the 
power spectrum. The power spectrum is found to be 
practically identical to that of ACDM for redshifts 
z > zssb, but as soon as the symmetry breaks at the 
background level, the symmetron fifth force can kick in 
and enhance the clustering of matter. 



It is clear from Fig. (12 1 that there exist a large 
range of parameters in which the symmetron model can 
be easily distinguished from ACDM. 

5.3. Halo profiles for x 



In Fig. ( 14 1 we show the profile of x inside the most 
massive halos found in the simulation. Since the fifth 
force is proportional to x this figure also provides infor- 
mation about the fifth force in halos. 

The field profile of x is seen to increase from the inner 
to the outer regions of the halos and shows that the fifth 
force is most suppressed in the central region as expected. 

The fifth force is stronger for smaller halos, because 
those generally reside in low density regions where the 
fifth force is less suppressed. We see the closer the sym- 
metry breaking redshift is to 0, the smaller x becomes 
inside the halo and the more suppressed the fifth force 
is. Again this is because early symmetry breaking means 
a higher critical density and halo needs to be more mas- 
sive to be effec tively screened. This effect is also seen 
3|9| 1 1 (notice the difference in distribution of 

"or 



in Fig. 

X at z 



between the differe nt ru ns) and also on the 
matter power spectrum in Fig. ( 12 1 (notice the way the 



power spectrum starts growing again on small scales for 
zssb = 2.0). 

This has some important consequ ence s for the local 
constraints. We mentioned in Sec. (2.1) that the local 
constraints were derived by assuming that our galactic 
neighborhood was not screened today, and lead to the 
constraint 



L(l + z SSB ) 3 / 2 <2.3, 



(61) 



From our numerical results we see that when zssb = 2.0 
its only the most massive halos that are screened. This 
means that the assumption that went into the constraint 
above is very likely to be true. On the other hand for SSB 
that happens very close to today, halos of much smaller 
mass are in fact screened and it might be possible to have 
a range L that exceeds Eq. (61 1 and still be in agreement 
with experiments. 

We note that we have not seen any significant effect 
on the halo density profiles. For a given mass range, 
the halo profiles seem to have approximately the same 
distribution. There should be some important differences 
for low mass halos, but the resolution in our simulation 
is too low to study this. 
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Fig. 7. — The fifth force to gravity for in a slice of the simulation box at z = (left) and the comparison between four different redshift 
(right) for a runs with zssb = 0.5, /3 = 1 and L = 1. Black shows 2 = 0, red z = 0.25, blue z = 0.65 and g reen z = 1. The solid lines shows 

F F ' t> - = 2/3 2 xi( z ) which is the prediction for the (short-range) forces at the perturbation level (see Eq. (pE)!). Note that the force will be 

dispersed around this prediction because we already have significant over and under densities in which Xr Xb- At z = 1 the background 
field is close to \b = and the force is small everywhere in space. As we move closer to z = 0, the symmetry breaks, and the background 
value moves towards \b = 1- This means that the force in low density regions (small gravitational force) will increase whereas in high 
density regions (strong gravitational force) the screening kicks in and the force becomes suppressed just as seen above. The numerical size 

H 2 

of the forces are given in terms of code units which are times the physical force unit. 



However, the halo number counts was significantly dif- 
ferent as we shall see in the next section. 

5.4. Halo mass function 

The halo mass function n is another key structure for- 
mation observable. It is defined to be the number density 
of dark matter halos within a given mass range. Because 
of the symmetron fifth force we expect more halos to be 
formed relative to the standard ACDM scenario. 

We first look at the total number of halos (the inte- 
grated mass function) with more than 100 particles which 
clearly shows the effect of the fifth force, see Table, pj. 

In Fig. (15 1 we have shown the mass function of the 
symmetron compared to ACDM at z = 0. We see a 
significantly higher mass function, especially for low mass 
halos which are generally found in low density regions 
where the fifth force is unscreened. The earlier symmetry 
breaking occurs and the stronger the coupling strength 
/3, the more halos are formed in agreement to what we 
would naively expect. 

The mass function converges to that of ACDM at very 
large halo masses for most parameters we have looked 
at. This is because the most massive halos have taken 
a very long time to form and therefore when the sym- 
metron kicks in at some low redshift the halo is already 
massive enough to be screened. However, for the largest 
z ssb — 2.0 we do have small increases in both the halo 
number density and the mass of the most massive halos. 
There have been reports of some tension between ob- 
servations and ACDM predictions with regards to very 
massive halos. Unfortunately for the symmetron model 
to be able to elevate this tension significantly we would 
need values of the parameters which are in conflict with 
local experiments. 



There is a large range of viable parameters for the sym- 
metron where the mass function deviates significantly 
from ACDM. 

6. SUMMARY AND CONCLUSIONS 

The symmetron mechanism is one way a scalar field 
can be non-minimally coupled to matter and still evade 
local gravity experiments. The symmetron model suffers 
the same fine tuning problems as chameleon models, but 
has the advantage of looking like a more natural effective 
theory. 

The energy density of the symmetron is too low to con- 
tribute to the dark energy and we must therefore add a 
cosmological constant to get accelerated expansion of the 
Universe. The background evolution of the symmetron 
model is simply indistinguishable from that of the ACDM 
model. 

This degeneracy is broken by the linear perturbations. 
In particular we have shown that the linear growth index 
7(2, k) can have a significant scale and redshift depen- 
dence together with a value today which can be distin- 
guished from the ACDM prediction for a large part of 
the parameter space. 

The structure formation in the nonlinear regime was 
investigated by using N-body simulations. N-body simu- 
lations have the advantage over linear theory in its abil- 
ity of fully capturing the nonlinear environmental de- 
pendence of the symmetron field. Our results confirm 
the expectation that in high density environments the 
fifth force becomes screened. Consequently, the key ob- 
servables such as the nonlinear matter power spectrum 
is closer to the ACDM predictions that expected from a 
linear analysis. 

We found that the symmetron can still produce large 
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Fig. 8. — The density distribution (left) and scalar-field distribution (right) for a run with zssb = 0.5, ft = 1 and L 
bottom z = 1, z = 0.66 and z = 0. 



1. From top to 



A 


B 


C 


D 


E 


F 


G 


H 


ACDM 


1634 


1694 


1678 


1871 


1758 


2051 


1671 


1788 


1607 



TABLE 2 

The halo count for our nine simulations. The corresponding symmetron parameters can be found in Table. |lj). 



observable signatures in both the nonlinear matter power 
spectrum and the halo mass function, which could in 
principle be detected by current and near future cosmo- 
logical observations such as Euclid. 

Note that in the simulations performed in this work, 
we have treated baryons as dark matter. However, since 
the symmetron field has a uniform coupling to all mat- 
ter fields we expect that all the results will qualitatively 
remain even after baryons are included. 

In conclusion, the symmetron model has been found to 
have a wide range of observable cosmological effects on 
both linear and nonlinear scales. This adds to the list 



of observational sig natures like making galaxies brighter 
(Davis et al. 20111 and the possibility o f being detected 
in near future local gravity experiments ( Hinterbichler & 
Khoury 2010 1 to mention some. The symmetron is there- 
fore a good candidate for the detection of new physics 
beyond the standard model. 
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APPENDIX A: USEFUL EXPRESSIONS 

Up to first order in the perturbed metric variables H, 
\1/ the non-zero components of the symmetric Lcvi-Civita 



connection are 



1 00 
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1 Ok 


= 2, k 


pi 
1 00 




pi 
1 Ofc 




r o 
L jk 




pi 
L jk 





(1 - 22- 2*) - 
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(66) 
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Fig. 10. — The density distribution (left) and scalar-field distribution (right) for a run with zssb = 2, /3 = 1 and L = 1. From top to 
bottom z = 2.33, z = 1, z = 0.66 and 2 = 0. 
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Fig. 11. — The full nonlinear power spectrum for a run with 
z S sb = L = j8 = 1.0 (dashed black) and ACDM (solid black). 
For comparison we also show the corresponding predications from 
linear perturbations theory in red. We clearly see the effectiveness 
of the screening mechanism. The linear predications do not 
take the symmetron mechanism into account and are hugely 
overestimating the power on small scales relative to ACDM. 

From these expression we find that the components of 
the Ricci tensor and Ricci scalar are found to be 



(68) 
(69) 
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tf" + -(S' + 3#') 



(70) 



(71) 



a* \ a 

APPENDIX B: DISCRETISATION OF EQUATIONS 

To implement the nonrelativistic equations into our nu- 
merical code, we have to rewrite then using code units, 
which are given by 



x 



H n B' 



tr 



(H B) 



Pc 



= tH , X 
V = BV X 



(72) 



where subscript c stands for code units, B is the boxsize, 
Ho = lOOkm/s/Mpc and an overline denotes background 
quantities. In what follows we shall write V = V c for 
simplicity. 

6.1. Scalar field equation of motion 
The equation of motion for \ m code units becomes 

™2 / x 2 



(BH y 



x + x 

2 



r) 



_p_ 

Ho 



M p i 

M 



(PcX - x) 



(73) 



where x is the background solutions and we have used 
<fio = x to simplify. Note that x varies in the region < 
X 2 < 1. Discreticed this equation becomes L h (xij,k) = 
where 



L> KXi,j,k) = Tl (r,u \1 U»+l,3",fc ~ 2 Xi,j,k + Xi-l,j,k) 
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'Mpi s " 
M 



3Q T 



{pcXi,j,k - x) 



(74) 



The Newton-Gauss-Seidel iteration says that we can ob- 
tain a new and more accurate solution of xfj"k usm S our 
knowledge about the old solution Xifk as 



xfTk = xt* k 



L h {xttk) 



dL h (xttk)/9xtf,k 



(75) 



where 



6 ac 



d Xi,, 



i,j,k 



h 2 (BH 



(76) 



so ( M ^ , 



6.2. Poisson Equation 

Since we can neglect the scalar field contribution to 
the Poisson equation, it remains unmodified from that of 
ACDM and reads (in code units) 



V 2 $ c = ~fl m (Pc,i,j,k ~ 1) 



6.3. Particle equation of motion 



(77) 



Using the code units, Eq. (56 1 can be rewritten as 



dXc _ Pc 

dt c a 2 
dt c 
The factor 



X T7* f (3M\ 2 f c'Wx d X \ 



(78) 



(79) 



can be also rewritten in terms of L, 



(3 and z S SB by using Eq. ( 22 1 
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Fig. 12. — The fractional difference in the nonlinear power spectrum relative to ACDM for {zssb = 1.0, L = 1.0} (top left), {zssb = 
1.0, £ = 2.0} (top right), {zssb = 0.5, £ = 1.0} (bottom left) and {zssb = 2.0, L = 1.0} (bottom right). For each case we show the results 
for the two values /3 = 0.5 and f3 = 1.0. 
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Fig. 13. — The fractional difference in the nonlinear power spectrum relative to ACDM for {zssb = 0.5, L = 1.0, (9 = 1.0} (left) and 
{ 2 SSB = 1.0, £ = 1.0, P = 1.0} (right) for several different redshifts. 
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